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Abstract 

A stochastic solution is constructed for a fractional generalization 
of the KPP (Kolmogorov, Petrovskii, Piskunov) equation. The solu- 
tion uses a fractional generalization of the branching exponential pro- 
cess and propagation processes which are spectral integrals of Levy 
processes. 

1 Introduction: The notion of stochastic so- 
lution 

The solutions of linear elliptic and parabolic equations, both with Cauchy and 
Dirichlet boundary conditions, have a probabilistic interpretation. These are 
classical results which may be traced back to the work of Courant, Friedrichs 
and Lewy [1] in the 1920's and became a standard tool in potential theory[2] 
[3]. For example, for the heat equation 

1 d 2 

d t u(t,x) = -— u(t,x) with u(0,x) = f(x) (1) 
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the solution may be written either as 



or as 

u(t,x)=EJ(X t ) (3) 
K x meaning the expectation value, starting from x, of the process 

dX t = dW t 

W t being the Wiener process. 

Eq.(l) is a specification of a problem whereas (2) and (3) are solutions in 
the sense that they both provide algorithmic means to construct a function 
satisfying the specification. An important condition for (2) and (3) to be 
considered as solutions is the fact that the algorithmic tools are independent 
of the particular solution, in the first case an integration procedure and in 
the second the simulation of a solution-independent process. This should 
be contrasted with stochastic processes constructed from a given particular 
solution, as has been done for example for the Boltzman equation[4]. 

In contrast with the linear problems, for nonlinear partial differential 
equations, explicit solutions in terms of elementary functions or integrals 
are only known in very particular cases. However, if a solution-independent 
stochastic process is found that (for arbitrary initial conditions) generates 
the solution in the sense of Eq.(3), a stochastic solution is obtained. In 
this way the set of equations for which exact solutions are known would be 
considerably extended. 

The stochastic representations recently constructed for the Navier-Stokes[5] 
[6] [7] [8] and the Vlasov-Poisson equations [9] [10] define solution-independent 
processes for which the mean values of some functionals are solutions to these 
equations. Therefore, they are exact stochastic solutions. 

In the stochastic solutions one deals with a process that starts from the 
point where the solution is to be found, a functional being then computed 
along the whole sample path or until it reaches a boundary. In all cases 
one needs to average over many independent sample paths to obtain a ex- 
pectation value of the functional. The localized and parallelizable nature of 
the solution construction is clear. Provided some differentiability conditions 
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are satisfied, the process also handles equally well simple or very complex 
boundary conditions. 

Stochastic solutions also provide an intuitive characterization of the phys- 
ical phenomena, relating nonlinear interactions with cascading processes. By 
the study of exit times from a domain they also sometimes provide access to 
quantities that cannot be obtained by perturbative methods[ll]. 

One way to construct stochastic solutions is based on a probabilistic in- 
terpretation of the Picard series. The differential equation is written as an 
integral equation which is rearranged in a such a way that the coefficients 
of the successive terms in the Picard iteration obey a normalization condi- 
tion. The Picard iteration is then interpreted as an evolution and branching 
process, the stochastic solution being equivalent to importance sampling of 
the normalized Picard series. This method is used in this paper to obtain a 
stochastic solution of a nonlinear partial differential equation, which is a frac- 
tional version of the Kolmogorov-Petrovskii-Piskunov (KPP) equation[12]. 



2 A fractional nonlinear partial differential 
equation 

We consider the following equation 

t Dy (t, x) = X - x D p e u (t, x) + u 2 (t, x)-u (t, x) (4) 

We use the same notations as in the study of the linear problem in [13]. t D" 
is a Caputo derivative of order a 

( 1 ft f (m) (r)d,T _ i . 

t D° /(*) = < r t~ a) j ° (*-^) Q+1 ~ m m i < a <m ^ 



d_ 



:f(t) a = m 



m integer. x Dg is a Riesz-Feller derivative defined through its Fourier symbol 
by 

T [ x Dgf (x) } (k) = (k) T {/ (x)} (k) (6) 

with^(fc) = |A;|' 3 e< si g n/£ W 2 . 

Eq. (4) is a fractional version of the KPP equation, studied by probabilistic 
means by McKean[14]. Physically it describes a nonlinear diffusion with 
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growing mass and in our fractional generalization it would represent the 
same phenomenon taking into account memory effects in time and long range 
correlations in space. 

As outlined in the introduction, the first step towards a probabilistic 
formulation is the rewriting of Eq.(4) as an integral equation including the 
initial conditions. For this purpose we take the Fourier transform (JF) in 
space and the Laplace transform (£) in time obtaining 

s a u{s,k) (s,k) = s a - 1 u(0 + ,k) — ^(k)u(s,k)-u(s,k)+ / dte~ st T (u 2 (t, x)) 

(7) 



where 



/oo 
e ikx u (t, x) dx 
-oo 

POO 

u(s,x) — C(u(t,x)) — / e~ st u(t,x)dt 

Jo 



This equation holds for < a < 1 or for < a < 2 with ^m(0 + ,x) = 0. 
Solving for u(s, k) one obtains an integral equation 

u(s,k) = : — a , - u (0 + , k) + / dt 1 — Q - - T (u 2 (t, x)) 

1 ; s a + 1 + m (k) K 1 J s a + 1 + (k) K 1 >> 

(8) 

Taking the inverse Fourier and Laplace [15] transforms 



u 



(...) = { -n £ ^ ta) ) (x - y) u (,+,,) 

+ f dr{t-T) a - l E^ a {-{t-T) a ) 
Jo 



f 

j — c 



^ r ( "l + (T^" rn '<'-^-'> (9> 



-E7a )P is the generalized Mittag-Leffler function 

oo _.j 



r (aj + p) 
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We define the following propagation kernel 
and, from the normalization relation, 



E aA (-t a )+ / drit-rT^E^i-it 



tT) 



o 



we may interpret E a ^ (—t a ) and (t — r) Q_1 i? Q)Q; (— (t — r) a ), respectively as 
a survival probability up to time t and as the probability density for the 
branching at time r in a branching process B a . It is a fractional general- 
ization of an exponential process. This provides a probabilistic sampling of 
the Picard series obtained by iteration of Eq.(9). The solution is therefore 
obtained by the expectation of the exit values of the following process: 

Starting at time zero, a particle lives according to the process B a . At the 
branching time r the initial particle dies and two new particles are born at 
the dying point. The process continues in the same way with independent 
evolution of each one of the newborn particles. At time t the solution is 
obtained as a functional of the n existing particles at time t, namely as the 
product of the initial condition propagated from the point where each one of 
the n particles is at time t up to the initial position. 



u(t,x)=E x (<p 1 <p 2 ---<p n ) (11) 



with 



<Pi = J dyfdyf ■ ■ ■ dy^_MKa (n,* - vf) G^ a (r 2 , y? - yf) ■ ■ ■ 
■ ■ ■ Gi, a (n-u yt 2 ~ vt.) <i (r k , ytl - yf) u (V, yf) (12) 

with Xlj=i Tj = t, k — 1 being the number of branchings leading to particle 
i. Notice that the last propagator in (12) is different from the others. 

Because of the normalization of the probabilities in the process B a , the 
probability of each one of the products in (11) corresponds to the weight of 
the corresponding term in the Picard series. Therefore the expectation value 
exists whenever the Picard series converges. 
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The solution (11) is not yet a purely stochastic solution because it involves 
both the expectation value over the process B a and a multiple integration of 
the initial condition with the propagation kernels G^ a l and G^ a . To obtain a 
purely stochastic solution we notice that, for < a < 1, the propagation ker- 
nels satisfy the conditions to be the Green's functions of stochastic processes 
in H. (see the Appendix). 

We denote the processes associated to G a l (t,x) and G@ (t,x), respec- 
tively by n^ ;1 and 11^ . Therefore the process leading to the solution is 
as described before with all the particles until the last branching propagat- 
ing according to the process and the last ones (that sample the initial 
condition) propagating by the process IT^. When finally all the n surviv- 
ing particles reach time zero, their coordinates x + £j are recorded and the 
solution is given by 

u(t, x) = E x (u(0 + , x + 6M0 + , x + 6) • • • u(0 + , x + Q) (13) 

Eq.(13) is a stochastic solution of (4) and our main result is summarized as 
follows: 

Theorem: The nonlinear fractional partial differential equation (4), with 
< a < 1, has a stochastic solution given by (13), the coordinates x + £j 
in the arguments of the initial condition obtained from the exit values of a 
propagation and branching process, the branching being ruled by the process 
B a and the propagation by 11^ 1 for the particles that reach time t and by 
Ha,a f or a H remaining ones. 

A sufficient condition for the existence of the solution is 

\u(0 + ,x)\<l (14) 

Remarks: 

1) The condition |w(0 + , x)\ < 1 imposes a finite value for all contributions 
to the multiplicative functional. However, the solution may exist under more 
general conditions, namely when the decreasing value of the probability of 
higher order products in (13) compensates the growth of the powers of the 
initial condition. 

2) The stochastic solution may also be constructed by a backwards-in- 
time stochastic process from time t to time zero. This is obtained by rewriting 
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Eq.(9) as 



r i (E a i{- (1 + Wb (*0) t a )\ 
u(t,x) = E a>1 (-t a ) J ^ dyF- 1 ( Q;1 V t (-t°) )( x -y) u (° + » 

+ / drr^E^ (-r a ) 
Jo 

7 oo ^ ( jL(-^) J (X ^ y) (t - T ' ^ (15) 

and noticing that also 

E a>1 (-t a ) + / drr^E^ (-r a ) = 1 
Jo 

Then, we obtain the following stochastic construction of the solution: 

Starting at time t a particle propagates backwards in time according to 
the process IT^ 1 if it reaches time zero or according to IT^ a if it branches at 
time t — r. The branching probability is controlled by the process B a (that 
is, the branching probability density is r a-1 i? a)a (— r a )). When it branches, 
two new particles are born which propagate independently and the process 
is repeated until all surviving particles reach time zero. 

Appendix. The Green's functions and the 
characterization of the processes 

The processes and H"^ a 

(*,*)} (t,k) = a ' 1{ { Eai{ 2 _/ a) y >} ' (16) 

r a , E a a (- (l + \tlA (k)) t a ) 

F{Gi, a M} (t,k) = M E aa (-£) (17) 

For a propagation kernel G (t, x) to be the Green's function of a stochastic 
process, the following conditions should be satisfied: 

(i) G(0,x-y) =5(x-y) or F {G} (0, k) = 1 VA: 

(ii) / dxG (t, x) = 1 Vt or T {G} (t, 0) = 1 
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(iii) G (t, x) should be real and > 
For the processes LT^ 1 and IT^ a 

(i) F (0, k) = J ! and T {GU (o, k ) = = 1 

W ^ (*. °) = = 1 and T i G U (*> °) = ferSr = 1 

(iii) If JF {G} (t, — fc) = {T {G} (t, k))* then G (t, x) is real. 
Because ^ (— fc) = (V^ (A;)) it follows 



and 



( - ( i + ^? (-*)) t a ) = (- (i + \tft (k) ) r 



implying that both G^ a l (t, x) and G^ a (t, x) are real. 

Finally, for the positivity, one notices that for < a < 1 and p > a, 
E a , P (~ x ) is a completely monotone function [16]. Therefore 



p oo 

E a , p {-x) = / e- rx dF(r) 
Jo 



with F nondecreasing and bounded. 

For G^ p (t, x) (p = 1 and p — a) one has 

GP(t,x) = r / dF(r) / dke- ihx e- rta ( 1+1 ^-V) 

= / dF{r)e- rta / dke- ikx e -?rW-Q 

We recognize the last integral (in k) as the Green's function of a Levy pro- 
cess. Therefore one has an integral in r of positive quantities implying that 
(*a,i x ) an d G^ a (t, x) are positive. 
The process B a 

The decaying probability in time dr of this process is 
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From ^ 

[ T«- x E a , a (-r a ) dr = l- E a>1 (-t a ) 
Jo 

it follows that E a> i (—t a ) is the survival probability up to time t. The process 
B a is a fractional generalization of the exponential process. 
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